


 % Related parameters
aux  = (1-alpha)/((sigma^2)/2);
saux = aux*s*lambda;
 

 % Solve for unknown constant of integration that first emerges in 
 % equ. for density of surplus in natural wastage region. 
 % Constants in other regions linked to this unknown. 
 % Integrate density for array of possible values 
 % for this unknown. Solution is value that
 % ensure density integrates to 1.
 %
 % To come up with a grid of constants, use
 % distribution function from natural 
 % wastage region
%G_mM  = (A / saux) *(mL^aux2) *((mM/mL)^aux2 -1). 
%G_mM  < 1 & G_mM >0 -->
Ahi  = .99*saux*(mL^-saux)/((mM/mL)^saux -1);
Alo  = .01*saux*(mL^-saux)/((mM/mL)^saux -1);
obs_const= 500; 
constvec = linspace(Alo, Ahi, obs_const)';
 %
 % 
 % Elsby/Gottfries show that unknown constant shows up
 % in BC b/c it enters the density of marg. product @ m=mL.
 % For any lambda and constant, U share of searchers
 % can be set to ensure BC is satisfied:
ushr  = (1/saux)*constvec*(mL^saux) ./ (1 + (1/saux)*constvec*(mL^saux));
 %
 %
 % Natural wastage (mL,mM)
obs_nw = 100;
[grid_nw, wghts_nw] = legendre(mL, mM, obs_nw);
G_mM = (constvec / saux) *(mL^saux) *((mM/mL)^saux -1);	
 % J(m) in NW region:
J_nw = -(omega0/rho0) + ((1-omega1)/rho1)*grid_nw + J1*(grid_nw.^gamma1) + J2*(grid_nw.^gamma2);
 %     
 %
 %
 % Full replacement (mM,mR)
obs_fr= 750; 
[grid_fr, wghts_fr] = legendre(mM, mR, obs_fr);
G_mR  = G_mM + (mM^saux)*constvec*log(mR/mM);
 % J(m) in FR region:
obs = obs_fr;
R_m = zeros(obs_fr, 1);
for j = 1:obs_fr
    m = grid_fr(j);
    [grid,wghts] = legendre(mM, m, obs);   
    delta    = s*lambda ./ (1 + aux*s*lambda*log(grid/mM));
    integrand=  (  scriptw *((m./grid).^phi1) + ...
                (1-scriptw)*((m./grid).^phi2) ).*(delta./grid);
    R_m(j) = c*aux * wghts'*integrand;	% replacement costs
end
J_fr = -(omega0/Rho0) + ((1-omega1)/Rho1)*grid_fr - R_m + JJ1*(grid_fr.^phi1) + JJ2*(grid_fr.^phi2); 
 %
 %
 % Net Expansion (mR,mU)
obs_ne  = 50;   % distance between mR and mU is very small but certain objects, such as eta, display *a lot* of curvature here.
[grid_ne, wghts_ne] = legendre(mR, mU, obs_ne);
 %
part1 = ((1-omega1)/(alpha*c))*(grid_ne - mR);
part2 = (r+((omega0+r*C)/c))*(log(grid_ne) - log(mR));
delta_mR = s*lambda ./ (1 + aux*s*lambda*log(mR/mM));
auxconst = ( ((1-omega1)/(alpha*c))*mR - (r+delta_mR + ((omega0+r*C)/c)) )*(mR^(-1/(1-alpha)));
part3 = auxconst*(1-alpha)*(grid_ne.^(1/(1-alpha)) - mR^(1/(1-alpha)));
 %
intdeltaoverm = (part1 - part2 - part3)*ones(1,obs_const);
G_ne   = (ones(obs_ne,1)*G_mR').*exp(aux*intdeltaoverm) + ...
         (ones(obs_ne,1)*(ushr./(1-ushr))').*(exp(aux*intdeltaoverm) - 1);  %obs_ne by obs_const 
 %
gconst = G_mR + (ushr./(1-ushr));      
gridmat_ne  = grid_ne*ones(1,obs_const);
deltamat_ne = ((1-omega1)/(alpha*c))*gridmat_ne - (r+((omega0+r*C)/c)) ...
             - auxconst*(gridmat_ne.^(1/(1-alpha)));
g_ne   = (ones(obs_ne,1)*gconst') .* aux.*(deltamat_ne./gridmat_ne) .* exp(aux*intdeltaoverm);     
 %
 %
 %
 % Sum up over three regions:
integral_of_g = G_ne(end,:)';   % Evaluate at mU.
if min(integral_of_g) >1 || max(integral_of_g) < 1
    disp('ERROR: density cannot integrate to 1');
    pause;
end
 % Solution for unknown constant and recover ushr (unemployed share of
 % searchers). With ushr, we can also compute unemployment and 
 % labor force.
constsol = interp1(integral_of_g, constvec, 1);
ushr     = (1/saux)*constsol*(mL^saux) ./ (1 + (1/saux)*constsol*(mL^saux));
urate    = s*ushr/(1 - ushr*(1-s));
L = En_firm/(1-urate);
 %
 %
 % Final density
g_nw = constsol*(grid_nw.^(saux-1));
g_fr = (mM^saux) * constsol*(grid_fr.^-1);  
g_ne = interp1(constvec, g_ne', constsol);	
g_ne = g_ne';	
if min(g_ne) <1e-9
    disp('ERROR: Density fnc. in NE region is negative');
    pause;
end
 %
 % Final c.d.f.
G_mM  = interp1(constvec,G_mM,constsol);
G_mR  = interp1(constvec,G_mR,constsol);
G_nw  = cumsum( wghts_nw.*g_nw );           % G_nw(end) does == G_mM
G_fr  = G_mM + cumsum( wghts_fr.*g_fr);     % G_fr(end) does == G_mR
G_ne  = G_mR + cumsum( wghts_ne.*g_ne);
 % 
wghts_all= [wghts_nw;  wghts_fr;  wghts_ne];
grid_all = [grid_nw; grid_fr; grid_ne];
g_all = [g_nw; g_fr; g_ne];
G_all = [G_nw; G_fr; G_ne];

